Graphic Solution of a Differential Equation of the First Order

Diffeq[x1_,y1_,f_,h_,n_]:=Module[{w=Table[0,{i,1,5},{j,1,n}]},w[[1,1]]=x1;w[[2,1]]=y1;Do[w[[3,i]]=f[w[[1,i]],w[[2,i]]];w[[1,i+1]]=w[[1,i]]+h;w[[2,i+1]]=w[[2,i]]+w[[3,i]]h,{i,1,n-1}]
;w];
f0[x_,y_]:=2y;
f1[x_,y_]:=y Log[y]^2;
f2[x_,y_]:=y/(x^2+1);

f3[x_,y_]:=3x y/(1+y);
f4[x_,y_]:=Sqrt[y/x];
f5[x_,y_]:=(2x y+y^2+-x^2+2)/2;
f6[x_,y_]:=y^2+x;
fs={f0,f1,f2,f3,f4,f5,f6};

Manipulate[With[{f=fs[[pf]]},With[{w=Diffeq[x1,y1,f,h,n]},With[{gr=Grid[{Join[{" "},Range[n]],Join[{Style["x",Italic]},w[[1,All]]],Join[{Style["y",Italic]},w[[2,All]]],Join[{Style["y",Italic]'},w[[3,All]]]},Frame->All]},
Text@Column[{TraditionalForm[y'==f[x,y]],Show[If[sn,With[{s=Quiet[NDSolve[{y'[x]==f[x,y[x]],y[x1]==y1},y,{x,x1,x1+(n-1)h}]]},Plot[Evaluate[{y[x]}/.s],{x,x1,x1+h(n-1)},PlotStyle->{Thickness[0.005],Brown},PlotRange->{{-.2,1.33},{-0.2,1.31}}]],Graphics[{}]],Graphics[{PointSize[0.013],{RGBColor[0, 0.866667, 0.866667],Thickness[0.008],Line[Table[{w[[1,i]],w[[2,i]]},{i,1,n}]]},If[shexp,{Red,Line[{{w[[1,n-1]],w[[2,n-1]]},{w[[1,n]],w[[2,n-1]]}}],Arrowheads[Small],Arrow[{{w[[1,n]],w[[2,n-1]]},{w[[1,n]],w[[2,n]]}}],Text[Style["h",Italic,10],{(w[[1,n-1]]+w[[1,n]])/2,w[[2,n-1]]},{0,1}],Text[Style[Row[{Style["h",Italic]," ",TraditionalForm[y'[(Subscript[x, n-1])]]}],10],{w[[1,n]],(w[[2,n-1]]+w[[2,n]])/2},2{-1,0}]},{}],Table[Text[Style[(Subscript[Style["x",Italic], i]),10],{w[[1,i]],0},3{0,1}],{i,1,n}],
Table[Text[Style[(Subscript[Style["y",Italic], i]),10],{0,w[[2,i]]},4{1,0}],{i,1,n}],Point[Table[{w[[1,i]],w[[2,i]]},{i,1,n}]],Table[Tooltip[Point[{w[[1,i]],w[[2,i]]}],
Style[Row[{(Subscript["P", i]),"(",(Subscript[Style["x",Italic], i]),", ",(Subscript[Style["y",Italic], i]),")"}],10]],{i,1,n}]}],AspectRatio->Automatic,ImageSize->{400,400},PlotRange->{{-.2,1.33},{-0.2,1.31}},Axes->True,GridLines->{Table[x1+i h,{i,0,n-1}],Table[w[[2,i+1]],{i,0,n-1}]},GridLinesStyle->GrayLevel[0.9],AxesOrigin->{0,0},ImagePadding->30],gr},Alignment->Center]
]]],

{{pf,1,Style["f",Italic]},Range[7]},
"",
"initial values",
{{x1,0.1,TraditionalForm[(Subscript[x, 1])]},0.1,0.2,0.01,Appearance->"Labeled", ImageSize-> Tiny},{{y1,.3,TraditionalForm[(Subscript[y, 1])]},.1,.6,0.01,Appearance->"Labeled", ImageSize-> Tiny},
"",
"step",
{{h,0.2,Style["h",Italic]},0.1,0.2,0.01,Appearance->"Labeled", ImageSize-> Tiny},
"",
"number of points",
{{n,3,Style["n",Italic]},2,6,1,Appearance->"Labeled", ImageSize-> Tiny},
"",
"show 'exact' solution",{{sn,False,""},{False,True}},
"",
"show explanation",{{shexp,False,""},{False,True}},
SaveDefinitions->True]

This Demonstration presents Euler's method for the approximate (or graphics) solution of a first-order differential equation with initial condition, y'=f(x,y), y((Subscript[x, 1]))=(Subscript[y, 1]).

The method consists of calculating the approximation of y(x) by:

(Subscript[x, i+1])=(Subscript[x, i])+h,

(Subscript[y, i+1])=(Subscript[y, i])+h y'((Subscript[x, i]))=(Subscript[y, i])+h f((Subscript[x, i]),(Subscript[y, i])),

where i=1,…,n-1.

These coordinates determine points (Subscript[P, 1]), (Subscript[P, 2]), …, (Subscript[P, n]). These points form Euler's polygonal line that is an approximate solution of the problem. We compare it with a better solution provided by Mathematica's built-in NDSolve function (brown line).

L. Euler, Institutionum calculi integralis volumen primum, section 2, chapter 7: De integratione aequationum differentialium per approximationem, written 1763, published 1768.

Contributed by: Izidor Hafner